Published online 17 January 2014 



Nucleic Acids Research, 2014, Vol. 42, No. 6 4031-4042 

doi:10.1093/nar/gktl388 



NAR Breakthrough Article 

A dynamic alternative splicing program regulates 
gene expression during terminal erythropoiesis 

Harold Pimentel 1 , Marilyn Parra 2 , Sherry Gee 2 , Dana Ghanem 2 , Xiuli An 3 , Jie Li 3 , 
Narla Mohandas 3 , Lior Pachter 1 ' 4 ' 5 and John G. Conboy 2 * 

department of Computer Science, University of California, Berkeley, CA 94720, USA, 2 Life Sciences Division, 
Lawrence Berkeley National Laboratory, Berkeley, CA 94720, USA, 3 Red Cell Physiology Laboratory, New York 
Blood Center, New York, NY 10065, USA, department of Mathematics, University of California, Berkeley, CA 
94720, USA and department of Molecular & Cell Biology, University of California, Berkeley, CA 94720, USA 

Received October 29, 2013; Revised December 16, 2013; Accepted December 17, 2013 



ABSTRACT 

Alternative pre-messenger RNA splicing remodels 
the human transcriptome in a spatiotemporal 
manner during normal development and differenti- 
ation. Here we explored the landscape of transcript 
diversity in the erythroid lineage by RNA-seq 
analysis of five highly purified populations of 
morphologically distinct human erythroblasts, rep- 
resenting the last four cell divisions before enucle- 
ation. In this unique differentiation system, we found 
evidence of an extensive and dynamic alternative 
splicing program encompassing genes with many 
diverse functions. Alternative splicing was particu- 
larly enriched in genes controlling cell cycle, 
organelle organization, chromatin function and 
RNA processing. Many alternative exons exhibited 
differentiation-associated switches in splicing effi- 
ciency, mostly in late-stage polychromatophilic 
and orthochromatophilic erythroblasts, in concert 
with extensive cellular remodeling that precedes 
enucleation. A subset of alternative splicing 
switches introduces premature translation termin- 
ation codons into selected transcripts in a differen- 
tiation stage-specific manner, supporting the 
hypothesis that alternative splicing-coupled 
nonsense-mediated decay contributes to regulation 
of erythroid-expressed genes as a novel part of the 
overall differentiation program. We conclude that a 
highly dynamic alternative splicing program in 
terminally differentiating erythroblasts plays a 
major role in regulating gene expression to ensure 
synthesis of appropriate proteome at each stage as 
the cells remodel in preparation for production of 
mature red cells. 



INTRODUCTION 

Alternative pre-messenger RNA (mRNA) splicing enables 
individual genes to generate multiple protein products that 
differ structurally and functionally by insertion or deletion 
of important functional domains encoded by alternative 
exons. During normal development and differentiation, 
dynamic changes in the expression or activity of the 
splicing regulatory machinery coordinately modulate 
networks of alternative splicing events in a spatiotemporal 
manner. Post-transcriptional RNA processing can thereby 
modulate essential protein functions according to the 
physiological requirements of the cell by regulating 
coherent biological pathways (1). Conversely, network 
perturbations caused by mutations that alter splicing 
factor expression or activity underlie an array of 
complex human diseases (2). The experimental work sup- 
porting these concepts has been performed primarily 
in non-hematologic cell types. However, recent studies 
have revealed that normal T cells execute a complex 
splicing program (3) and that splicing factor mutations 
are associated with hematological cancers, including 
myelodysplasia (4-6), demonstrating the importance of 
alternative splicing in hematology. Characterization of 
the alternative splicing program in human erythroblasts 
undergoing terminal erythroid differentiation will reveal 
novel insights into RNA regulatory pathways that drive 
cell differentiation and provide a basis for identifying 
splicing defect in human erythroid diseases. 

Alternative isoforms of various erythroid transcripts 
have been reported, and a modest number of differenti- 
ation stage-specific switches in alternative splicing patterns 
are known (7). Best studied mechanistically is the upregu- 
lation of splicing efficiency for protein 4.1R alternative 
exon 16, which encodes part of the spectrin-actin 
binding domain required for optimal assembly of a mech- 
anically competent membrane skeleton (8,9), an essential 
structure for mature erythrocyte function. This exon is 
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predominantly skipped in early erythroblasts but included 
efficiently in late erythropoiesis (10,11), and several RNA 
binding proteins/splicing factors that influence exon 16 
splicing have been identified (12-15). Underscoring the 
importance of this splicing switch, failure to include 
exon 16 causes mechanically unstable red cells and 
aberrant elliptocytic phenotype with anemia (16,17). 

To comprehensively explore the alternative splicing 
landscape during terminal erythroid differentiation, 
we used an RNA-seq strategy to analyze and compare tran- 
scriptomes of highly purified human erythroblasts cultured 
in vitro from CD34 + cord blood progenitors (18). RNA-seq 
allows a robust high-resolution assessment of the transcrip- 
tome, but there remain computational challenges in data 
interpretation. Using extensions of current transcript abun- 
dance estimation tools combined with non-parametric stat- 
istical methods, we found an extensive alternative splicing 
program that is significantly enriched in genes controlling 
cell cycle, organelle organization, chromatin function and 
RNA processing. Importantly, hundreds of these alterna- 
tive splicing events are regulated in a differentiation stage- 
specific manner, with most switches in exon inclusion/ex- 
clusion efficiency occurring in late-stage erythroblasts con- 
current with extensive cellular remodeling as erythroblasts 
transition from a highly proliferative state to a terminally 
differentiated state. Finally, we discovered a subset of 
splicing switches that introduce premature translation ter- 
mination codons (PTCs), thus decreasing the proportion of 
full-length coding mRNAs and downregulating gene ex- 
pression via induction of nonsense-mediated decay 
(NMD). Post-transcriptional RNA processing pathways 
may be regulated by this mechanism in late erythroid dif- 
ferentiation, as numerous RNA binding proteins exhibit 
elevated expression of PTC exons in the most mature 
erythroblasts. These studies support the central hypothesis 
that human erythroblasts execute a surprisingly complex, 
differentiation stage-specific alternative splicing pro- 
gram that is essential for normal proliferation and 
differentiation. 

MATERIALS AND METHODS 

Human erythroblast populations 

CD34+ cells were enriched from cord blood, cultured under 
conditions that promote erythroblast differentiation, 
and further enriched for stage-specific erythroblast popula- 
tions by fluorescence-activated cell sorting (FACS) exactly 
as described (18). By gating narrow windows, high cell 
purity is attained for five discrete erythroblast populations 
(>90%), representing the last four cell divisions in terminal 
erythropoiesis. In some experiments, the cell populations 
enriched for proerythroblasts (proE) (culture day 9) and 
orthochromatic erythroblasts (orthoEs) (culture day 16) 
were treated with lOOug/ml cycloheximide for 4h to 
inhibit nonsense-mediated decay of RNA transcripts con- 
taining PTCs. 

RNA-seq analysis 

Total RNA was extracted using an RNeasy Plus Mini Kit 
(QIAGEN) from sorted human erythroblasts at distinct 



stages of erythroid maturation. The integrity of the 
RNA was evaluated on an Agilent Bioanalyzer, and 
samples with RNA integrity numbers (RIN) >8.0 were 
prepared for sequencing. Poly(A)-tailed RNA was subse- 
quently prepared by the Epigenomics Core, Cornell 
Medical College, using the mRNA Seq Sample Prep Kit 
(Illumina Inc., San Diego, CA, USA) and used to create 
libraries for HiSeq2000 sequencing (Illumina). 

Hundreds of millions of 50-nt RNA-seq reads were 
obtained from each of the 15 erythroblast samples, repre- 
senting five distinct stages of differentiation from three 
independent sorting experiments. RNA-seq fragments 
were aligned to the Ensembl-annotated transcriptome 
version 67 (19) [hgl9 human genome build GRCh37 
(20)] using the Bowtie aligner (21). Read statistics are 
given in Supplementary Table SI. Transcript-level esti- 
mates were obtained using the transcript abundance 
estimation tool, eXpress (22), and were expressed as 
fragments per kilobase of exon per million fragments 
mapped (FPKM) as described previously (23,24). For 
each exon of every Ensembl gene, transcripts from the 
gene were divided into two sets: the set of isoforms 
including the exon (termed inclusion) and the set 
of isoforms with an intronic region spanning it (termed 
overlap). Expression of individual exons in 'exon- 
inclusion' isoforms relative to the total expression of all 
isoforms was represented as percent spliced in (PSI) and 
computed by calculating the ratio FPKM (inclusion)/ 
FPKM (inclusion, overlap). Because only ratios of reads 
were compared, a library size normalization step between 
samples was not needed. 

To classify exons as differentially expressed between the 
proE and orthoE samples, we used the estimated PSI 
values with the samr R package (25), assuming continuous 
data. This resulted in a test similar to a Mest, but where 
the null distribution was generated by a permutation test 
non-parametrically from the data. In the same manner, we 
estimated the PSI of each exon in the basophilic (basoE) 
and polychromatophilic (polyE) erythroblast populations. 
We then considered any exon for which the inclusion plus 
overlap FPKM was in the lower quartile to be insuffi- 
ciently expressed and discarded it from the analysis. We 
also removed any exon for which a contributing isoform 
had a contradiction in coverage by having long spanning 
regions of 0 coverage (see Supplementary Methods). From 
this list filtered for low gene expression, we also filtered 
out exons with minimal alternative splicing (PSI <0.05 or 
>0.95) at any differentiation stage, because empirically we 
found that this PSI filter greatly reduced the incidence 
of false positives. Data for the PSI-filtered set of splicing 
events are shown in Supplementary Table S2, and 
the larger unfiltered PSI data set is provided in 
Supplementary Table S3. The filtered list of alternatively 
spliced exons was used as input into GOrilla {Gene 
Ontology enRIchment anaLysis and visuaLizAtion tool) 
for Gene Ontology (GO) analysis (26) (Table 1) and for 
analyzing stage-specificity of splicing switches (Figure 5). 

Besides the computationally defined alternative splicing 
switches, a few additional examples were discovered in 
genes of interest by manual inspection of RNA-seq 
reads mapped to the human genome. Exons that did not 



Nucleic Acids Research, 2014, Vol. 42, No. 6 4033 



Table 1. GO terms enriched in alternatively spliced erythroblast genes 
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computationally meet a APSI threshold of 10% but 
showed strong evidence for stage-specific changes in the 
ratio of exon junction reads supporting inclusion versus 
exon exclusion reads were tested by reverse-transcriptase 
polymerase chain reaction (RT-PCR); examples of 
splicing switches detected in this manner are MEF2D, 
P1CALM and MTMR3. 



Validation of computationally predicted alternative exons 

Many alternative splicing events were verified by manual 
inspection of mapped reads on the human genome 
browser, to filter out exons not supported by multiple 
exon-inclusion and exon-skipping reads. In selected 
cases, RT-PCR analysis was performed to visualize the 
amplified products corresponding to inclusion and 
skipping events. Total RNA was reverse-transcribed into 
cDNA using random primers and the Superscript III First 
Strand Synthesis System (Invitrogen), and PCR analysis 
was performed using AccuTaq polymerase (Sigma) using 
primers located in constitutive exons located upstream 
and downstream of each candidate alternative exon. 
Primer sequences are given in Supplementary Table S4. 
For initial RT-PCR validations, RNA from unsorted 
day 9 (proE-enriched) and day 16 (orthoE-enriched) cells 
was used; these results were subsequently confirmed with 
RNA purified from highly enriched (FACS-sorted) proEs 
and orthochromatophilic erythroblasts. 



RESULTS 

Differentiating erythroid cells execute a robust alternative 
splicing program 

Highly purified human erythroblasts were obtained from 
in vitro cultures by FACS using differentiation stage- 
restricted cell surface markers as previously described 
(18). Stage-specific RNA-seq libraries were prepared 
from five morphologically distinct populations of succes- 
sively more differentiated erythroblasts, each separated 
from the next by one cell division: proEs, early basophilic 
erythroblasts (e-basoE), late basophilic erythroblasts 
(1-basoE), polychromatophilic erythroblasts (polyE) and 
finally orthochromatophilic erythroblasts (orthoE), the 
last stage before enucleation (Figure 1). Three independ- 
ent preparations of sorted cells were analyzed to provide 
biological replicates. 

Alternative splicing events in differentiating erythro- 
blasts were identified following the strategy shown in 
Figure 2. A total of 1.5 billion reads representing the tran- 
scriptomes of all 1 5 erythroblast samples were mapped to 
the Ensembl annotated transcriptome and analyzed using 
numerous tools to derive exon-level alternative splicing 
information as described in 'Methods' section. The 
initial set of 243 464 candidate alternative exons was 
filtered to remove low-expression genes, and exons with 
minimal alternative splicing, to derive a data set of 3784 
alternative exons in 2250 genes (Supplementary Table S2). 
Validation of computationally predicted alternative 
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Figure 1. Morphology of differentiation stage-specific erythroblast populations purified by FACS as described previously (18). 
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Figure 2. Flow chart showing the process by which major alternative 
splicing events were detected in discrete populations of erythroblasts at 
specific stages of their terminal differentiation. 
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Figure 3. RT-PCR validation of alternative splicing events in human 
erythroblasts. Shown are amplification products containing (filled 
circles) or lacking (empty circles) alternative exons amplified from 
erythroblast transcripts. Gene names are indicated above each lane. 



splicing events was performed by manual inspection of 
RNA-seq reads mapped to the human genome and 
by RT-PCR analysis of the relevant transcript regions. 
Figure 3 illustrates representative alternative splicing 
events in genes involved in transcriptional regulation 
(NFE2L1, MEF2D, SOX6, ATRX, FOXM1, MAX), 
DNA and histone modifications (DNMT1, MBD1, 
CARM1, KDM4B) and signal transduction pathways 
(MAP2K7, MAP3K3, MAP3K7). The majority of these 
alternative exons did not substantially alter their splicing 
efficiency between proEs and orthoEs; i.e. similar results 
were observed by computational and experimental ana- 
lyses independent of differentiation stage. However, alter- 
native exons in NFE2L1, MEF2D and CARM1 were 
among those observed to splice differently as the cells 
matured (discussed below in Figure 4). Of interest, 
SOX6 and ATRX are important regulators of globin 
gene expression (27-29). These results demonstrate that 
differentiating erythroblasts execute an extensive alterna- 
tive splicing program that can regulate protein structure 
and function via control of pre-mRNA splicing; 
moreover, when the affected proteins are transcription 
factors, important secondary effects on transcript levels 



of other genes might amplify downstream cellular 
consequences. 

To investigate more broadly which biological processes 
are impacted by the erythroid splicing program, we per- 
formed GO analysis. GO terms significantly enriched 
among genes expressing major alternative splicing 
events, relative to all genes expressed (above the lowest 
quartile) in erythroblasts, are shown in Table 1. The top 
terms included numerous descriptors associated with chro- 
matin structure and function, cell cycle regulation, organ- 
elle organization and regulation of RNA splicing. Among 
the read-validated alternative splicing events, we found at 
least 24 differentially expressed exons in 22 genes encoding 
enzymes that participate in histone modifications (30) 
(Supplementary Table S5), and at least 65 high-confidence 
alternative splicing events were detected in 48 genes 
encoding RNA binding proteins (Supplementary 
Table S6). As indicated in the tables, many of these alter- 
native splicing events are not annotated among Refseq 
transcripts, and ~40% are predicted to introduce PTCs 
that induce NMD. All of these splicing events passed low- 
expression and low-alternative splicing filters, and all were 
validated by RNA-seq read mapping. These observations 
suggest a high level of complexity in regulation of 



Nucleic Acids Research, 2014, Vol. 42, No. 6 4035 



A 

proE 

Differentiation-associated 
PSI increase 

orthoE 



gene: EPB41 DMTN NDEL1 NEK1 CARM1 BNIP2 NFE2L1 
stage: P 0 P 0 P 0 P O P 0 P 0 P O 




4 



PSI: 20 59 30 64 6 52 28 70 30 52 28 70 27 66 

gene: MEF2D* US01 MTMR3* MARK2 LGALS9 HPS1 PICALM* AXIN1 
stage: p O P Q P Q P Q P O P O P O P O 




PSI: 19 26 5 47 41 39 21 55 14 61 28 44 50 59 28 70 



B 



Differentiation-associated 
PSI decrease 

P'° E f I f 1 

orthoE 1— [] ^- 

gene: HSF2 
stage: P O 




PSI: 69 23 



C Differentiation-associated switch 
of mutually exclusive exons 

p'° E — A^Ph 

orthoE 1— |p[ f \ 

H2AFY 
P O 



PSI: 63 22 
25 60 



Figure 4. Differentiation stage-specific alternative splicing switches in human erythroblasts. (A) Alternative exons that exhibit stage-specific 
upregulation in splicing efficiency. Gels show splicing changes assessed by RT-PCR using primers located in flanking constitutive exons. 
Amplification products from proE RNA and orthoE RNA are shown in the left and right lanes, respectively. Gene name is indicated above 
each gel, while the calculated PSI is shown below each lane. Arrowheads indicate exon inclusion products. Asterisk indicates genes for which 
PCR results indicate larger splicing changes than were predicted bioinformatically. (B) Alternative splicing event that exhibits a differentiation- 
associated decrease in PSI as cell matures. (C) A differentiation-associated switch involving mutually exclusive exons of 100 and 91 nt, respectively. 
The size difference allows electrophoretic separation of the alternative products. 



erythroid gene expression, and a strong potential that 
these splicing events have downstream effects on the abun- 
dance and structure of transcripts of numerous other 
genes. Thus, important biological pathways express a 
rich transcript diversity during late erythropoiesis. 

Stage-specific switches in alternative splicing remodel the 
transcriptome during late erythroid differentiation 

Alternative splicing switches can play key roles during 
differentiation and development by selective up- or 
downregulation of individual exons that alter function of 
the encoded protein(s). To determine the extent of splicing 
transitions during erythroid differentiation, we first 
compared PSI values for alternative exons in early 
(proE) and late (orthoE) erythroblast populations. In all, 
1166 exons for which the difference in PSI values (APSI) 



exceeded 10% were identified by this approach; 266 of 
these exhibited APSI values exceeding 25%. A subset of 
splicing switches that alter protein coding domains was 
validated by RT-PCR using primers in the flanking con- 
stitutive exons (Figure 4). The previously described switch 
in protein 4.1R (EPB41) pre-mRNA splicing served as a 
positive control; as expected, splicing of exon 16 exhibited 
a substantial increase in PSI in orthoE cells compared with 
a much lower splicing efficiency in proE (Figure 4A). 
Similar to 4.1R, many of the newly discovered switches 
involve upregulation of splicing in orthoE (Figure 4A), 
but examples of downregulation (Figure 4B) and 
switching between mutually exclusive exons (Figure 4C) 
were also observed. For some of these, e.g. NDEL1, 
USOl, MEF2D and MTMR3, the 'new' isoforms with 
potentially new functions can be generated in orthoE 
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Table 2. Alternative splicing switches that alter coding domains 
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that were absent or poorly expressed in proE. Except for 
USOl, PICALM 1 and LGALS9, most of the genes with 
upregulated coding exons were also upregulated at the 
gene expression level, relative to global gene expression 
levels, further supporting the idea that the new isoforms 
might play important functions in late erythroblasts. 

Consistent with a role for alternative splicing in import- 
ant erythroid pathways, several validated switches in 
alternative splicing map to genes that function in cytoskel- 
etal organization, cell division and chromatin function 
(Table 2). One of the biggest switches occurred in the 
NDEL1 (nuclear distribution factor E-homolog-likel) 
gene that is important for nuclear migration and 
nucleokinesis (31), while a less dramatic switch was pre- 
dicted in survivin (BIRC5), a cytokinesis factor with a 
novel role in erythroblast enucleation (32,33). In the 
H2AFY gene, relative expression of mutually exclusive 
100- and 91-nt exons reversed almost completely 
between proE and orthoE (Figure 4A). This switch 
alters the structure and function of the modified histone, 
macroH2Al, which can replace the canonical H2A to gen- 
erate functional differences in chromatin activity (34-36). 
Consistent with the H2AFY splicing switch in terminal 
erythroid differentiation, differences in the relative abun- 
dance of the two isoforms have been reported in other 
tissues rich in proliferating versus differentiating cells 
(37). Other switches of likely functional importance 
occurred in kinases having roles in chromosome segrega- 
tion and cytokinesis [NEK1, (38^10)] and microtubule- 
and actin-based cytoskeletal networks [MARK2, (41)] 
in a clathrin assembly gene (PICALM) associated with 
anemia and iron metabolism abnormalities in mouse 
mutants (42,43); and the mitotic assembly protein 
NUMA1. 

To explore the hypothesis that alternative splicing 
networks can impact erythroblast remodeling, we 
examined the differentiation stage-specificity of splicing 



transitions. PSI as a function of erythroblast stage is 
shown in Figure 5A for several PCR-validated alternative 
exons that exhibited substantially increased inclusion 
between proE and orthoE (APSI values >34). For these 
genes, most of the change in alternative splicing occurred 
in the last two (most mature) erythroblast populations, 
although lesser transitions in splicing efficiency occurred 
at earlier stages. Extending the analysis to the broader set 
of alternative exons, we next determined APSI values 
across each of the four-stage transitions: proE to early 
basoE, early to late basoE, late basoE to polyE and 
polyE to orthoE. Figure 5B shows that ~155 stage- 
specific splicing switches with APSI > 25 were identified 
in these cultures, and the great majority of these 
occurred during the last two cell divisions. Of note, this 
is a conservative estimate that counts only exons having 
PSI values between 5 and 95 in all 15 samples assayed. 
These results demonstrate that most splicing transitions 
occur in late erythroblasts, temporally correlated with ex- 
tensive cellular remodeling as the cells prepare to enucle- 
ate, and supporting the hypothesis that splicing changes 
play a determinative role in promoting late erythroid 
differentiation. 

Alternative splicing coupled to NMD regulates expression 
of selected RNA processing factors 

Alternative splicing coupled to NMD (AS-NMD) can 
regulate gene expression by splicing-mediated introduc- 
tion of PTCs that induce nonsense-mediated degradation 
of the affected transcripts. PTC-containing transcripts 
can be generated by inclusion of exons that encode stop 
codons or alter translational reading frame (PTC-upon- 
exon inclusion), or by exon skipping events that alter 
reading frame (PTC-upon-exon skipping). Having noted 
that a number of RNA binding proteins express PTC tran- 
scripts in erythroblasts (Supplementary Table S4), we won- 
dered whether AS-NMD was functioning to maintain 
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homeostatic control of expression as described in previous 
studies (44^6). Instead, PSI and gene expression data 
indicated a different situation characterized by substantial 
stage-specific changes in relative expression of PTC 
transcripts (Figure 6). Increased proportions of PTC- 
containing transcript isoforms were expressed in late 
erythroblasts by either increases in exon inclusion or by 
increases in exon skipping (Figure 6A). Many of the pre- 
dicted increases are substantial; e.g. the SNRNP70 gene 
switches from PSI = 13 to PSI = 43 for the smaller PTC- 
upon-inclusion exon. Similarly, well-known AS-NMD 
targets in the SR protein family including SRSF2, 
SRSF3, SRSF6 and SRSF7 also showed large APSI 
switches from proE to orthoE (APSI values of +12 to 
+47). Other classes of AS-NMD candidates include 
hnRNP proteins (D, L, LL and M), small nuclear ribo- 
nucleic particle proteins (SNRPA1, SNRNP70 and 
U2AF1), alternative splicing regulators [TRA2A (45), 
TRA2B] and a 3' processing factor (CPSF4). In a few 
cases, an increased proportion of PTC transcripts in late 
erythroblasts can occur in the context of PTC-upon-exon 
skipping events (Figure 6A, lower panel). SNRPA1, for 
example, expressed mainly full-length transcripts in proE 
as indicated by PSI = 94, but the decrease to PSI = 53 in 
orthoE predicts substantial expression of PTC isoforms in 
orthoE. However, it is important to note that some known 



[SRSF11, SNRPB, (45)] or candidate (TIAL1 and 
HNRNPD) PTC splicing events showed no evidence of 
stage-specific regulation (Figure 6B, upper panel), and a 
few (RHD, BCL2L12, CLK1 and CLK4) were regulated 
so as to decrease the proportion of PTC isoforms in 
orthoE (Figure 6B, lower panel). RT-PCR validation of 
splicing switches involving PTC exons enriched in orthoE 
or proE is provided in Figure 7. We propose that splicing- 
mediated generation of non-productive (PTC-containing) 
transcripts is an important post-transcriptional regulatory 
mechanism for controlling a selected subset of RNA 
processing factors in late-stage erythroblasts. 

In principle, elevated expression of PTC transcripts 
could be due to increased production (splicing regulation) 
or to decreased degradation (NMD regulation). We 
reasoned that simple loss of NMD activity could not 
explain maturation-associated increases in PTC exons, as 
only selected PTC exons were upregulated in orthoE and a 
few were actually downregulated in these cells. This 
hypothesis was investigated experimentally by assaying 
whether inhibition of NMD could affect steady-state 
levels of PTC exons. If NMD is active in normal 
orthoE, then steady-state PSI values for PTC exons 
should increase due to selective stabilization of these tran- 
scripts in cycloheximide-treated erythroblasts inhibited for 
NMD. Figure 8 shows that PTC-upon-inclusion exons in 
SNRNP70 and HNRPLL had substantially elevated PSI 
values in NMD-inhibited cells, with PSI of a 72-nt exon in 
SNRNP70 increasing moderately in proE-enriched day 9 
erythroblasts and dramatically in orthoE-enriched day 16 
erythroblasts (upper panel). NMD inhibition also 
increased the proportion of PTC transcripts from the 
BCL1L12 and CLK4 genes that were generated by exon 
skipping and were more enriched in proEs (Figure 8, lower 
panel). We conclude that the NMD machinery is still 
active in late erythroblasts and therefore that the accumu- 
lation of PTC isoforms in these normal cells is not primar- 
ily due to loss of NMD, but instead must be due, at least 
in part, to regulation at the level of alternative splicing. 
Moreover, as the 'missing' quantity of PTC isoforms in 
normal cells must have been degraded, we propose that 
alternative splicing coupled to NMD (AS-NMD) is a 
mechanism for downregulating gene expression in late 
erythroblasts, analogous to the intron retention mechan- 
ism used by late-stage granulocytes to downregulate 
selected genes (47). Consistent with this, 9 of 11 RBP 
genes with upregulated PTC events (82%) were 
downregulated at the gene expression level relative to 
total mRNAs in late erythroblasts. In contrast, among 
genes shown in Figure 4 with validated upregulation of 
non-PTC exons in orthoE, 12 of 15 genes (80%) exhibited 
increases in gene expression level. 



DISCUSSION 

The erythroblast system provides a unique opportunity to 
study the role of post-transcriptional RNA processing in a 
differentiation stage-specific manner during the extensive 
cellular remodeling that characterizes terminal differenti- 
ation. We found a robust alternative splicing program 
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enriched in genes that function in cell cycle regulation, 
organelle organization, chromatin structure and function 
and RNA processing, all of which undergo tremendous 
changes before enucleation. Most importantly, dynamic 
switches in alternative exon usage remodel the late eryth- 
roid transcriptome so as to generate new mRNA isoforms 
or to alter proportions of existing isoforms. These data 
support the hypothesis that a highly structured alternative 
splicing program is executed during terminal erythroid 
differentiation to regulate key biological processes as the 
cells prepare to enucleate and form mature red cells. 

Previously, proper alternative splicing of protein 4.1 
(EPB41) pre-mRNA was shown to be essential for 
generating a red cell membrane with the specialized mech- 
anical properties required for survival in the circulation 
(48). The finding of abundant alternative splicing events 
in erythroblasts suggests many additional contributions to 
erythroid-specific functions. For example, regulators of 
alpha- and beta-globin gene expression such as ATRX 
and its interacting partners MECP2 and macroH2Al, as 
well as BCL11A and SOX6 all express multiple transcript 
isoforms whose importance has not been sufficiently 
evaluated. The GO categories enriched in alternative 
splicing events suggest additional roles in regulation of 



chromatin condensation, autophagy and enucleation in 
late erythroblasts, all key and unique events during late 
stages of terminal erythroid differentiation. We anticipate 
that some of these alternative splicing switches might be 
shared by other cell types undergoing the transition from 
proliferation to terminal differentiation; e.g. differences in 
relative abundance of H2AFY isoforms in tissues rich in 
proliferating versus differentiating cells (37) mirror the 
changes observed in erythroblasts. The erythroid differen- 
tiation system may be ideal for studying such splicing 
events. 

Our studies suggest that a second major function for the 
erythroid alternative splicing program is post-transcrip- 
tional downregulation of gene expression in erythroblasts. 
We hypothesize that NMD-associated splicing switches, 
similar to switches in coding exons, are driven by 
changes in expression of splicing regulatory proteins 
during terminal erythroid differentiation that can up- or 
downregulate expression of PTC exons, according to 
physiological requirements for the encoded proteins. The 
marked changes in steady-state levels of PTC-containing 
transcripts during differentiation indicate that the mech- 
anism regulating their production must be distinct from 
homeostatic (often autoregulatory) control mechanisms 
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previously ascribed to AS-NMD [(49) and references 
therein]. By maintaining active machinery for both alter- 
native splicing and NMD, erythroblasts could coordinate 
these processes to limit expression in proliferating proE 
of differentiation functions required for remodeling in 
orthoE, while also downregulating in orthoE proliferative 
functions not needed in late erythroblasts. Differentiation- 
associated switches in PTC exons have already been 
proposed to operate in differentiating neurons, as, for 
example, changes in expression of splicing factor PTB 
can regulate switches in both PTC and non-PTC splicing 
events (50,51). The full impact of such pathways on eryth- 
roid biology may extend well beyond RNA binding 
proteins, as a substantial fraction of alternative splicing 
events affecting histone modifying enzymes (~42%, 
Supplementary Table SI) and many DNA binding 
proteins (results not shown) were predicted to induce 
NMD. Moreover, the current study likely underestimated 
or entirely missed PTC splicing switches that generate 
rapidly degraded transcripts, or splicing switches that 
resulted in retained introns; the latter occurs during late 
stages of granulocyte differentiation where they can 
suppress expression of selected genes (47). In fact, inspec- 
tion of RNA reads mapping to the SF3B1 gene re- 
vealed increasing abundance of an intron retention event 
in late-stage erythroblasts (results not shown) that predicts 
reduced expression of this critical splicing factor 
implicated in myelodysplasia (4,52). Future studies will 
explore these potential expansions of AS-NMD by 
RNA-seq analysis of NMD-inhibited erythroblasts. 

In conclusion, coordinated temporal control of alterna- 
tive splicing plays a critical role in regulating structure and 
function of the erythroid transcriptome, likely having a 
major impact on orderly differentiation of erythroblasts. 
Analogous cell type-specific or differentiation-associated 
alternative splicing programs that are critical for cell 
identity have been identified in non-erythroid cell types 
including neurons, muscle cells, epithelial cells and 
T lymphocytes (3,53-55). Conversely, abnormal splicing 
factor expression as a consequence of human genetic 
disease or as a result of genetic knockdown in model or- 
ganisms can profoundly alter splicing networks and have 
adverse phenotypic effects. Disturbance of programmed 
splicing networks underlies many complex human 
genetic diseases including myotonic dystrophy, the auto- 
immune neurologic disease paraneoplastic opsoclonus- 
myoclonus ataxia, facioscapulohumeral muscular 
dystrophy and amyotrophic lateral sclerosis (53,56-59). 
Moreover, a number of hematological cancers including 
myelodysplasia can be caused by mutations in splicing 
factor genes (4,52,60,61), and it has been shown that 
erythroblasts from myelodysplasia patients exhibit dis- 
ordered erythropoiesis (62). Understanding the erythro- 
blast splicing program will provide new insights into 
normal differentiation mechanisms and also into dis- 
ordered differentiation in disease states. 
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